PRO BGPS_FFORE_GRID, block, subset
  
  IF n_elements(subset) EQ 0 THEN subset = ''
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  NBLOCK = 32
  
  IF n_elements(block) EQ 0 THEN BEGIN
     message,'Enter block number to grid (1-'+string(NBLOCK,format="(I0)")+$
             '):',/inf
     modstr = ''
     read,modstr
     block = fix(modstr)
  ENDIF
  IF block LT 0 || block GT NBLOCK THEN $
     message,'Block # is out of range.  Please try again...'
  
  
  
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_MODEL2, BGPS_FFORE_SIMPLE, $
     BGPS_FFORE_SIMPLE2, BGPS_FFORE_LOOPIE, STAR_DENSITY
  
  rotc = 0b

  ;;=====================================================================
  ;; Print system time for timing purposes
  print,'Start Time:' & spawn,'date'
  
  ;; Read in the IDL save file RB3
  IF n_elements(fn) EQ 0 THEN fn = 'bgps_rb3'+subset+'.sav'
  IDIR = './irdc_dist_model/'
  message,'Restoring RB3 file '+fn,/inf
  restore,IDIR+fn,/ver
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  ;; Do some work, and place this info into the structure!
  n       = n_elements(rb3)
  tau     = dblarr(n)
  f_data  = dblarr(n)
  sig_f   = dblarr(n)
  
  ;; Do the re-usable calculations
  GENERATE_LPSTR, ROTC=rotc
  do_Tfit = 0b

  message,'Starting the gridding...',/inf
  
  ;;================*****===========
  ;; Enter present number of models (or many more) -- only generates
  ;;   an array of filenames.
  
  
  message,'Case V:  (veil + gas + stars), rho_star fixed',/inf
  savfn = IDIR+'data/case3_irdc_grid'+subset+'_block'+$
          string(block,format="(I0)")+'.sav'
  
  ;; A is sampled every 0.01
  n_A = 51L
  A = dindgen(n_A)/100.d
  
  ;; B is sampled logarithmically
  n_B = 41L
  B = 10.d^(dindgen(n_B)*0.1d - 2.d)
  
  ;; C is sampled logarithmically
  n_C = 41L
  C = 10.d^(dindgen(n_C)*0.1d - 2.d)
  
  ;; Divide Td into blocks
  ;; Td is sampled every 0.25 K from 10 K to 30 K
  n_Td = 81L
  Td = dindgen(n_Td)*0.25d + 10.d
  
  i0 = long(round(float(block - 1)/float(NBLOCK)*float(n_Td)))
  i1 = long(round(float(block)/float(NBLOCK)*float(n_Td)))-1
  print,i0,i1
  
  Td = Td[i0:i1]
  n_Td = n_elements(Td)
  
  
  
  ;; Rest of the stuff...
  grid = fltarr(n_A,n_B,n_C,n_Td)
  nfar = intarr(n_A,n_B,n_C,n_Td)
  total_steps = n_A * n_B * n_C * n_Td
  steps = 0L
  t_start = systime(1)
  
  FOR m=0L, n_Td-1 DO BEGIN
     TEMP_DEPENDENT_VARIABLES, Td[m]
     
     FOR i=0L, n_A-1 DO BEGIN
        FOR j=0L, n_B-1 DO BEGIN
           FOR k=0L, n_C-1 DO BEGIN
              
              ;; Do the calculation
              p = [A[i],B[j],C[k],1.,Td[m],2.,0.]
              resid = BGPS_FFORE_LOOPIE(p,NEARFAR=nearfar)
              grid[i,j,k,m] = TOTAL(resid*resid,/NAN)
              nfar[i,j,k,m] = fix(TOTAL(nearfar))
              
              ;; Print information
              steps++
              pstr = string([A[i],B[j],C[k],Td[m]],$
                            format="(F0.2,'  ',F0.3,'  ',F0.3,'  ',F0.1)")
              frac = double(steps)/double(total_steps)
              etd = t_start + (systime(1) - t_start)/frac
              print,pstr+'    ['+systime(0,etd)+']'
              
           ENDFOR
        ENDFOR
     ENDFOR
  ENDFOR
  
  save,grid,nfar,filename=savfn
  
END
